Symmetric three-particle motion in Stokes flow: equilibrium for heavy spheres in 

contrast to "end-of-world" for point forces 
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A stationary stable solution of the Stokes equations for three identical heavy solid spheres falling 
in a vertical plane is found. It has no analog in the point-particle approximation. Three spheres 
aligned horizontally at equal distances evolve towards the equilibrium relative configuration while 
the point particles collapse onto a single point in a finite time. 



In recent experiments, the structure and velocity field 
of a sedimenting non-Brownian suspension have been 
observed to differ significantly from the equilibrium 
state 0. Theoretical explanation of this behavior must 
take into account hydrodynamic interactions between the 
suspended particles. The simplest interesting cluster con- 
sists of three identical point particles settling under grav- 
ity in a quiescent viscous infinite fluid. Its evolution was 
analyzed numerically and turned out to be sensitive to 
small changes of the initial configuration 0. Symmet- 
ric periodic motions of three point-particles and of three 
spheres located at vertices of an isosceles triangle have 
been also found 0,0] and the equilateral horizontal trian- 
gle is known as the equilibrium solution. Until now, theo- 
retical results confirmed the analogy between the motion 
of three well-separated spheres and three point-particles. 

In this work a simple example is studied, which illus- 
trates that evolution of well-separated spheres may dif- 
fer significantly from the corresponding point-particle ap- 
proximation. Settling of three identical particles under 
gravity is analyzed. The particle centers are located at 
vertices of an isosceles vertical triangle with the horizon- 
tal base. First, the system of equations is specified, and 
the method to solve it is outlined. Then, a new equilib- 
rium relative configuration of the spheres is found. Next, 
the particles are initially aligned horizontally and the mo- 
tion of the sphere centers is determined and compared 
with the motion evaluated within the point-particle ap- 
proximation. Finally, other relative trajectories of the 
spheres are determined and stability of the equilibrium 
is shown. 



sphere centers ri(t) satisfy the following equations, 



A low-Reynolds-number fluid flow is considered. Its 
velocity v(r) and pressure p(r) satisfy the Stokes equa- 
tions 0, 



r?V 2 v 



V P = 0, 



(1) 



with the fluid viscosity r). The fluid is infinite. Its motion 
is generated by settling of three identical spheres under 
gravity F = —Fe z . The stick boundary conditions are 
assumed at the surfaces of the spheres. Positions of the 



Lk=l 



F, 



i = 1,2,3. 



(2) 



The 3x3 mobility matrices fj, ik depend on the instan- 
taneous relative configuration of all the spheres, and are 
evaluated numerically by the multipole expansion 0, Q . 
The algorithm from Ref. Q and its accurate numerical 
FORTRAN implementation described in Ref. [|| are ap- 
plied, with the multipole order L = 4. The set of the or- 
dinary differential equations is solved numerically by 
the adaptive fourth-order Runge-Kutta method [ljj. I n 
the following, distances will be normalized by the sphere 
diameter d and time by two Stokes times t 8 = 3irr]d 2 / F , 
keeping the same notation. The dimensionless variables 
satisfy Eq. @ with F = 1. 

Initially, the sphere centers are located at vertices of an 
isosceles vertical triangle with the horizontal base. From 
the symmetry of the Stokes equations it follows that the 
time-dependent configuration is also of this type. The 
apex sphere is labeled 2 and the two other spheres are 
labeled 1 and 3. The relative positions r\ 2 =T\—r-x and 
= r3 — r 2 are parameterized as 



ri2 = (-x/2,0,z), 
r 32 = (x/2,0,z). 



(3) 
(4) 



The distance x between the twin spheres 1 and 3 and 
the vertical separation z between the twins and the apex 
sphere 2 satisfy the following system of equations, 



x = v x (x,z), 
z = v z (x,z), 



(5) 
(G) 



with the initial values a;(0) = xq, z(0) = zq. Here v x = 

J2k=li^3k,xz — l*2k,xz], and V z = J2 k=1 [^3k,zz — M2fc,zz], 

with the Cartesian components of \x ik dependent only 
on x and z. Once x(t), z(t) are evaluated, then r 2 = 
[0, 0, z 2 {t)\ is obtained by a direct integration of Eq. J5J. 

In the point-force approximation, the same units are 
used, and a single point moves with the Stokes velocity 
of the single sphere. The point-particle dynamics reads, 



ik ' &z &z 



(7) 



2 



with the dimensionless Oseen tensor, 

T lk = JL(I + f ifc f tt ), (8) 

the unit vectors = rife/Vjfc and the unit tensor 
I. Eqs. Q) are integrated numerically by the adaptive 
fourth-order Runge-Kutta method [l(| . 

Now, equilibria for the dynamical system ©-© will 
be found. The r.h.s. of Eq. JSJ) vanishes if the twin 
spheres 1 and 3 touch each other and the distance be- 
tween their centers x = 1. In this case, the lubrication 
forces prevent spheres 1 and 3 from a relative mo- 
tion. Therefore v x (l,z) = at each point z, and the 
motion is vertical. At equlibrium point z, also v z (l,z) 
should vanish. Consider first configurations with the sin- 
gle sphere 2 below the touching dublet. In this case z is 
positive; moreover, z > V3/2, because the spheres do not 
overlap. The function v z (l,z) is evaluated numerically 
and plotted in Fig. ^ Two equilibria are seen. The first 
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FIG. 1: Motion of the touching spheres with respect to the 
single one. 

one, with (x, z) — (l,v / 3/2), corresponds to the touch- 
ing triplet of spheres, with relative motions excluded by 
the lubrication forces It is interesting to see that 

there exists also another positive root z eq of the equation 
v z {\, z) — 0. For a given multipole order L, z eq is evalu- 
ated numerically by the standard bisection method [Ic) . 
Next, z eq is evaluated for all L = 3, ...,30 and the limit 
L — > oo is taken, 

z eq = 1.578634. (9) 

The question arises does this equilibrium attract tra- 
jectories. First, vertical motion of the horizontal touch- 
ing doublet located above the singlet will be analyzed, 
using the function plotted in Fig. ^ In the limit of infi- 
nite z, the faster doublet and the slower singlet are not 
influenced by each other and as a result i = —0.3799554, 



see e.g. [12|. Starting from any z > z eq , the relative 
distance z between the doublet and the singlet decreases 
with time, because the vertical velocity v z is negative. 
However, if zq < z eq , then v z > and z increases with 
time. It is remarkable that in this case a heavier doublet 
is repelled by a lighter singlet located below. In both 
cases, the system tends to the equilbrium position z eq , 
approaching it after infinite time. Indeed, from Fig.^it 
is clear that dv z {l 1 z) / dz < 0, if z is close to z eq , and 
therefore | z — z eq | decreases with time exponentially. 

Consider now configurations where the touching 
spheres are below the single one, with z < 0. The Stokes 
equations are invariant with respect to the time reversal, 
supperposed with the reflection in the horizontal plane, 
which contains the center of sphere 2. Therefore 

v x (x,z) = -v x (x, -z), v z (x,z) = v z (x, -z). (10) 

In particular, for z < the relative velocity is immedi- 
ately obtained from Fig.Q] using the relation v z (l, z) — 
v z (l, —z). Starting from negative Zq, trajectories escape 
from — z eq . If zq > —z eq , then the doublet is attracted 
by the singlet and z — > — \/3/2. If zq < —z eq , then the 
doublet is repelled by the singlet and z — > — oo. To con- 
clude, the dynamics i = v z (l, z) has four equilibria, with 
x = 1 and z = ±z eq or z — ±\/3/2. 

Now, general case will be studied when spheres 1 and 
3 do not touch each other. Assume first that initially 
three sphere centers are aligned horizontally. In Fig. |3 
typical trajectories Zi(xi) of spheres and point-particles 
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FIG. 2: Trajectories of three sphere centers (dashed lines, 
circles) and three point-particles (solid lines, stars), at t = 
centered at xi = 2 = —£3, X2 — 0, 21 = 22 = 23 = 0. The 
symbols represent positions at times t=10, 20 and 30. 



3 



are compared. Here the solution of the dynamics J2J is 
denoted as ri(t) = [xi(t),0, Zi(t)], with i = 1,2,3. At 
the beginning, trajectories of the point particles coincide 
with those of the spheres. However, after a finite time 
t = 36.74 they collapse onto a single point, and obviously 
the point-force approximation breaks down. 

In the following, only the relative two-dimensional dy- 
namics ©-© will be discussed. Evolution of the initial 
values xq and zq — is of special interest. The distance 
£ = x — 1 between the surfaces of spheres 1 and 3 de- 
creases monotonically, decaying exponentially to zero for 
long times, as illustrated in Fig-El F° r close twin spheres, 




FIG. 3: Time-dependent size £ of the gap between surfaces of 
the spheres 1 and 3; initially £o = 3 and zo = 0. 

if Xo is sufficiently small, then their vertical separation z 
from sphere 2 increases with time. For a large xq, at the 
beginning z increases, significantly exceeding z eq , then 
drops down below z eq to approach it slowly again, as de- 
picted in Fig. 0] Numerical integration of Eqs. ©-© 
is limited to the gap sizes £ < 10~ 14 . It is remarkable 
that for such close distances, the vertical separation z 
still differs from its equilibrium value z eq by around 4%. 
This difference is related to rotation of spheres 1 and 2, 
significant even for such tiny gaps between their surfaces. 

Close to the equilibrium, for very small £ and small 
z — z eq , it is possible to solve the dynamics ©-(jGl), using 
the asymptotic lubrication expansion of the two-sphere 
mobilities for spheres 1 and 3. To extract the leading 
terms in the limit of £ — > 0, the velocity of each sphere 1 
and 3 is written as superposition of velocities calculated 
for two other problems in the absence of sphere 2: sedi- 
mentation of the spheres 1 and 3 and the free motion of 
the spheres 1 and 3 in the ambient flow determined by 
the dynamics JSJ)-©- On the other hand, the free mo- 
tion is obtained by applying the two-sphere mobility to 
the forces Fpi and torques Tpi exerted on spheres i = 1, 3 




FIG. 4: Evolution of the vertical separation z(t) from the twin 
spheres 1,3 to sphere 2 (numerical result). Initially xo = 4, 
zo = 0. Inset: comparison with the asymptotic expression. 



fixed in the same ambient flow. Therefore 



X\ 



(y?i + yi 3 )(Fou - 1) + Hvn + y? 3 )T iy, 



(ii) 

(12) 



with the components of the two-sphere mobility denoted 

as in Ref. [ill b Y x iu x i3> Viv Vi3> Vli and Via- The force 
and torque units are F and Fd, respectively. For sphere 
3, the corresponding formulas follow from the symmetry 
with respect to reflections in the plane x — 0. For £ — > 0, 
the forces F i = F 03 and torques T i = — T 3, as well as 
the velocity of sphere 2, are regular functions of £. From 
Ref. [HI it follows that (^-4) ~ £, while {ytr+Vvs) ~ 
l/(ln£ + £) + const and (yf x +yf 3 ) - l/(ln£ + S), with 
B — 4.00. Therefore for very small £, 



v x (x,z) 
v z (x,z) 



f(z)£ + 0(e InO, 
h(z) 



ln£- 



B 



O 



(13) 
(14) 



The leading terms in the expansion of z around 
z eq give f(z) ~ —A, g(z) w —C(z — z eq ) and 
h{z) « —D. The constants A = 0.80, C = 0.124 
and D = 0.22 are evaluated numerically as A = 
-(dv x (x,Zeq)/dx) x =i, C = -(dv z (l, z)/dz) z = z<sq and 
D = -{dv z (x, z eg )/i9[l/(ln^ 1 + b)]} x=1 . Asymptotic 
approximation of Eqs. JSJ-lJBJ reads 



e = -M, 

z = -C[z - z eq ) 



D 



lnir 1 + B' 

By integrating Eq. I|15fl and the trajectory equation, 

D 



dz_ C_ 
dt~ Ai {z ' 



Zeq) 



^(ln^" 1 +B) 



(15) 
(16) 

(17) 
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the asymptotic solution is easily obtained, 



{(z a -Z c 



(D/A)[Ei(r) 



(18) 
(19) 



where r = C(ln£ _1 + B)/A = C[t+ (ln^ 1 +B)/A], and 
its value at t = is equal to t = C(ln CcT 1 +-B) /A. As be- 
fore, £o = Xo — 1. The symbol Ei(r) denotes the exponen- 
tial integral, the Cauchy principal value of J_ dte t /t. 

Precision of the asymptotic solutions (|18[1 - (|19[1 is illus- 
trated in the inset of Fig. For gaps £ < 10~ 4 , the 
numerical z(t) is approximated with 5% accuracy by the 
asymptotic one. Similar estimation holds for ln£(t). 

For very large times (for very small gaps £), t — * oo, 
and Ei(r) ~ (1 + l/r)e r /r. In this case, relation 1|19|) 
further simplifies, 



D 



C ln£" 



1 



B 



if£-0, 



(20) 



and (x, z) — > (x eg , z e9 ), if t — > oo. The equilibrium solu- 
tion (x eqi z eq ) of the dynamics {S}-|(BJ is stable. 

Finally, in Fig. we numerically evalute the phase 
portrait for the dynamics 101-® • The two-dimensional 
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FIG. 5: Phase portrait of sphere dynamics JSJ-©. The stable 
node is denoted by •. Inset: sphere trajectory (solid line) is 
compared with its point-particle counterpart (dashed line). 

phase space excludes non-overlapping spheres, and is 
therefore given as {(x, z) : x > 1 and x 2 /4 + z 2 > 1}. 
The equilibrium (x eq ,z eq ) is a stable improper node. At 
this point, dz/dx = — oo for all the trajectories. Indeed, 
dz/dx ~ -{D/C)/[{hi^- x + B) 2 £] close to (x eqi z eq ), as 



it follows from Eq. H20|) . The trajectories of the spheres 
initially aligned horizontally keep the same shape as long 
as all the spheres are at large distances from each other. 
As illustrated in the inset of Fig. \E\ this "ideal" shape is 
given by the scale- free point-particle dynamics. In Fig. [SI 
another family of trajectories is also plotted, with z — oo 
at t = — oo. These trajectories never reach z = 0. 

For z < 0, trajectories and the motion are easily ob- 
tained by the time reversal superposed with the reflec- 
tion in the plane 2 = 0, using Eq. I|10l) . In particular, 



(Xeqi 



j) is unstable improper node. Summarizing, the 



equilibrium points of the sphere dynamics jSj-ljSJl are lo- 
cated on the curve z 2 +x 2 /4 = 1 (marked by a dotted line 
in Fig. |SJ) and at (x eq ,±z eq ). On the contrary, no equi- 
librium exists for the relative motion of point-particles, 
as it follows from Eq. Q. 

To conclude, a stable equilibrium has been found for a 
symmetric relative motion of three identical spheres set- 
tling under gravity, in which the sphere centers form an 
isosceles vertical triangle with the horizontal base. In the 
equilibrium configuration, the two upper spheres touch 
each other and are well-separated from the lower sin- 
glet. For a large class of initial conditions, including dis- 
tant spheres, the system evolves towards the equilibrium, 
reaching it after infinite time. The corresponding point- 
force approximation breaks down, because after a finite 
time, all the particles experience "the end-of-world" , col- 
lapsing onto a single point. If two spheres are far below 
the single one, the system separates, increasing the rela- 
tive vertical distance with time. In addition to the physi- 
cal importance of the results, the solutions presented here 
may be used as benchmarks for numerical simulations of 
many-particle systems. 
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